rng('default');


% simulate to get unconditional joint kappa/lastIssuance density
Binit           = 10000;  % number of initial periods to discart
Ninit           = 100;    % number of firms
Tinit           = 120000;
[~,states]      = sim_markov_chain_fast(P_nodis,Tinit+Binit,1,1,1:nS-1);
g_sim           = mu_X(states) + sig_XS(states).*randn(Tinit+Binit,1) + sig_Xid(states).*randn(Tinit+Binit,Ninit); % [Tsim,Nsim]
k_def_sim       = k_def(states); % [Tsim,1]
k_iss_sim       = k_iss(states); % [Tsim,1]
k_bar_sim       = k_bar(states); % [Tsim,1]
k_sim           = NaN(size(g_sim));
k_QML           = NaN(size(g_sim));
last_issuance   = ones(size(g_sim)); % times initial state of the Markov chain
k_sim(1,:)      = k_bar(states(1));
k_QML(1,:)      = k_bar(states(1));
for t = 1:Tinit+Binit-1 % maturity
    k_lag                = k_sim(t,:);
    G_lag                = exp(g_sim(t,:));
    default              = k_lag <= k_def_sim(t);
    issuance             = k_lag >= k_iss_sim(t);
    k_sim(t+1,:)         = k_lag.*G_lag;
    k_sim(t+1,default)   = k_bar_sim(t).*G_lag(default);
    k_sim(t+1,issuance)  = k_bar_sim(t).*G_lag(issuance);
    k_QML(t+1,:)         = k_QML(t,:);
    k_QML(t+1,issuance)  = k_bar_sim(t);
    last_issuance(t:end,issuance)  = ones(Tinit+Binit-t+1,sum(issuance)).*states(t); % most recent issuance state
end
k_pdf     = k_sim(Binit+1:end,:);
k_QML_pdf = k_QML(Binit+1:end,:);
LI_pdf    = last_issuance(Binit+1:end,:);
k_pdf     = k_pdf(:);
k_QML_pdf = k_QML_pdf(:);
LI_pdf    = LI_pdf(:);
n_pdf     = length(k_pdf);
clear Binit Ninit Tinit states g_sim k_def_sim k_iss_sim k_bar_sim last_issuance k_sim t k_QML


% generate crisis paths, evaluate firm stats
NOBS_before  = 1;  % in years
NOBS_after   = 20; % in years
Tsim         = (NOBS_before+NOBS_after)*12;

Psim         = 100; % # of paths per rep
Nsim         = 100; % # of firms
Rsim         = 10000; % # of reps

MOM0 = NaN(Rsim,Psim);
MOM1 = NaN(Rsim,Psim);
MOM2 = NaN(Rsim,Psim);
MOM3 = NaN(Rsim,Psim);
parfor rr=1:Rsim
    
    % simulate aggregate paths
    epsE_paths   = NaN(Tsim,Psim);
    state_paths  = NaN(Tsim,Psim);
    duration_rr  = NaN(Psim,1);
    end_D        = NaN(Psim,1);
    C_loss       = NaN(Psim,1);
    for i=1:Psim
        
        s0          = sum(rand>cumsum(stat(1:nS-2)))+1;
        [~,states1] = sim_markov_chain_fast(P,NOBS_before*12,1,s0,1:nS);
        [~,states2] = sim_markov_chain_fast(P,NOBS_after*12,1,3,1:nS);
        states      = [states1(:);states2(:)];
        epsCE       = randn(Tsim,2)*chol([1,rho_CE;rho_CE,1]);
        epsC_sim    = epsCE(:,1);
        epsE_sim    = epsCE(:,2);
        delC_sim    = mu_C(states) + sig_C(states).*epsC_sim;
        
        num_D            = find([0;diff(states2(:))],1)-1;
        if isempty(num_D)
            num_D = length(states2);
        end
        C_loss(i)        = sum(delC_sim(NOBS_before*12+(1:num_D)));
        duration_rr(i)   = num_D;
        state_paths(:,i) = states;
        epsE_paths(:,i)  = epsE_sim;
        
    end
    duration(:,rr) = duration_rr;
    
    
    % simulate kappa panels
    g_sim                = mu_X(state_paths) + sig_XS(state_paths).*epsE_paths + sig_Xid(state_paths).*randn(Tsim,Psim,Nsim);
    k_def_sim            = k_def(state_paths); % [Tsim,Psim]
    k_iss_sim            = k_iss(state_paths); % [Tsim,Psim]
    k_bar_sim            = k_bar(state_paths); % [Tsim,Psim]
    init                 = round(rand(Psim,Nsim)*n_pdf+0.5);
    default              = false(Tsim,Psim,Nsim);
    k_sim                = NaN(Tsim,Psim,Nsim);
    k_sim(1,:,:)         = k_pdf(init);
    k_QML                = NaN(Tsim,Psim,Nsim);
    k_QML(1,:,:)         = k_QML_pdf(init);
    last_issuance        = NaN(Tsim,Psim,Nsim);
    last_issuance(1,:,:) = LI_pdf(init);
    
    for t = 1:Tsim-1
        
        kt    = squeeze(k_sim(t,:,:)); % [Psim,Nsim]
        qmlt  = squeeze(k_QML(t,:,:)); % [Psim,Nsim]
        kbt   = k_bar_sim(t,:)'*ones(1,Nsim); % [Psim,Nsim]
        gt    = squeeze(exp(g_sim(t,:,:))); % [Psim,Nsim]
        
        kt1   = kt.*gt;
        qmlt1 = qmlt;
        kbt1  = kbt.*gt;
        
        default_t        = kt <= k_def_sim(t,:)';
        kt1(default_t)   = kbt1(default_t);
        qmlt(default_t)  = kbt(default_t);
        default(t,:,:)   = default_t;
        
        issuance_t       = kt >= k_iss_sim(t,:)';
        kt1(issuance_t)  = kbt1(issuance_t);
        qmlt(issuance_t) = kbt(issuance_t);
        
        k_sim(t+1,:,:)   = kt1;
        k_QML(t+1,:,:)   = qmlt;
        
        st                     = state_paths(t,:)'*ones(1,Nsim);
        lit                    = squeeze(last_issuance(t,:,:));
        lit(issuance_t)        = st(issuance_t);
        lit(default_t)         = st(default_t);
        last_issuance(t+1,:,:) = lit;
        
    end
    
    
    % consumption drop
    MOM0(rr,:) = C_loss;
    
    % leverage
    Fit     = griddedInterpolant({1:nS,k},lamDx);
    DX_sim  = Fit(repmat(state_paths,[1,1,Nsim]),k_QML);
    Fit     = griddedInterpolant({1:nS,k},lamSx);
    SX_sim  = Fit(repmat(state_paths,[1,1,Nsim]),k_sim);
    QMLEV_sim = DX_sim./( DX_sim + SX_sim .* k_sim./k_QML);
    QMLEV_sim(default) = NaN;
    avg_QMLEV  = nanmean(QMLEV_sim,3);
    
    Fit        = griddedInterpolant({1:nS,k},IV_ATM);
    IV_sim     = Fit(repmat(state_paths,[1,1,Nsim]),k_sim);
    avg_IV     = nanmean(IV_sim,3);
    
    Fit        = griddedInterpolant({1:nS,k,1:nS},CDS(:,:,:,60));
    CDS60_sim  = Fit(repmat(state_paths,[1,1,Nsim]),k_sim,last_issuance);
    avg_CDS60  = nanmean(CDS60_sim,3);
    
    for xx=1:Psim
        ix               = duration_rr+1:Tsim;
        avg_QMLEV(ix,xx) = NaN;
        avg_IV(ix,xx)    = NaN;
        avg_CDS60(ix,xx) = NaN;
    end
    
    MOM1(rr,:) = 100*max(avg_QMLEV);
    MOM2(rr,:) = 100*max(avg_IV);
    MOM3(rr,:) = 100*max(avg_CDS60);
    
end

MOM0 = reshape(MOM0,[Rsim*Psim,1]);
MOM1 = reshape(MOM1,[Rsim*Psim,1]);
MOM2 = reshape(MOM2,[Rsim*Psim,1]);
MOM3 = reshape(MOM3,[Rsim*Psim,1]);

cut = 0:-0.025:-0.25;
LEV = NaN(length(cut)-1,1);
VOL = NaN(length(cut)-1,1);
CDS = NaN(length(cut)-1,1);
for i=1:length(cut)-1
    ix     = (exp(MOM0)-1)<=cut(i) & (exp(MOM0)-1)>cut(i+1);
    LEV(i) = mean(MOM1(ix));
    VOL(i) = mean(MOM2(ix));
    CDS(i) = mean(MOM3(ix));
end

